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We propose a novel algorithm with a modified Tucker decomposition for tensor network that allows for 
efficiently and precisely calculating the ground state and thermodynamic properties of two-dimensional (2D) 
quantum spin lattice systems, and is coined as the fast linearized tensor renormalization group (fLTRG). Its 
amazing efficiency and precision are examined by studying the spin- 1/2 anisotropic Heisenberg antiferromag- 
net on a honeycomb lattice, and the results are found to be fairly in agreement with the quantum Monte Carlo 
calculations. It is also successfully applied to tackle a quasi-2D spin- 1/2 frustrated bilayer honeycomb Heisen- 
berg model, where a quantum phase transition from an ordered antiferromagnetic state to a gapless quantum 
spin liquid phase is found. The thermodynamic behaviors of this frustrated spin system are also explored. The 
present fLTRG algorithm could be readily extended to other quantum lattice systems. 
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Efficient and accurate numerical methods are very cru- 
cial to tackle the strongly correlated quantum lattice sys- 
tems. Although some analytical techniques and numerical 
methods have been proposed in the past decades, a large 
class of intriguing correlated electron and spin models are 
still intractable owing to the complexity of quantum many- 
body systems. Several numerical renormalization group 
(RG) approaches were thus developed, where the density 
matrix renormalization group [1| and its finite temperature 
variant — the transfer matrix renormalization group [if] have 
achieved a great success for one-dimensional (ID) systems. 
Very recently, generalizing the RG-based algorithms to two- 
dimensional (2D) quantum lattice systems has been remark- 
ably advancing. A few numerical approaches, for instance, 
the projected entangle pair state (PEPS) Jst], the tree tensor 
network (TN) |0L the multiscale entanglement renormaliza- 
tion ansatz state||5tL the infinite PEPS 130], the tensor renor- 
malization group Il8l 4l0ll . and so on, were proposed, some 
of which already gained interesting applications (e.g. Refs. 
IUJ,LL2|]). It is noted that most of these algorithms are effec- 



tive for the ground state properties, but they are still difficultly 
applied to study the thermodynamics of 2D quantum lattice 
models. 

By incorpor ating the infinite time-evolving block decima- 
tion technique 111311 . we developed a linearized TRG (LTRG) 
algorithm IU4I1 that renders a convenient way to investigate the 
thermodynamic properties of low-dimensional quantum spin 
lattice systems. Although the LTRG method is quite efficient 
and accurate for ID quantum systems, its cost is relatively 
high and the performance near a critical point needs careful 
improvements for 2D quantum systems. Within the frame- 
work of LTRG, when the density operator is represented by 
a TN through Trotter-Suzuki decomposition [15], the trunca- 
tion is needed to prevent from the divergence of dimension of 
Hilbert space during the imaginary time evolution, which will 
unavoidably bring errors that become worse in 2D quantum 
systems. To solve this issue, we note that the Tucker decom- 
position (TD) 11161 Il7ll is a nice way to obtain the best lower 



dimensional approximation of a single tensor, and has wide 
applications in areas of data compression, image processing, 
etc. IU7I1 The algorithms for the TD like the higher-order sin- 
gular value decomposition 1 18 , 3 and higher-order orthogo- 
nal iteration (HOOI) lEoll were suggested. 

In this work, by extending the HOOI scheme to a TN in- 
stead of a single tensor for an optimal truncation, we propose 
a novel algorithm that allows us to efficiently and accurately 
simulate not only the ground state but also thermodynamic 
properties of 2D quantum spin lattice systems in the thermo- 
dynamic limit, which is dubbed as the fast LTRG (fLTRG). 
We find that the computational cost of fLTRG is insensitive to 
the coordination number without losing the accuracy, which 
allows for a higher bond dimension cutoff D c when it is ap- 
plied to 2D and quasi-2D quantum systems. The cost of the 
fLTRG algorithm is ~ 0(zD 3 c ), while the LTRG is ~ 0(23^ -3 ), 
where z is the coordination number. The reliability, efficiency 
and accuracy of the fLTRG algorithm are examined by study- 
ing a spin- 1/2 anisotropic Heisenberg antiferromagnet on the 
honeycomb lattice, whose energy per site, staggered magneti- 
zation and specific heat are efficiently and accurately obtained 
by the fLTRG, and the results are in good agreement with 
quantum Monte Carlo (QMC) results. To show its powerful 
performance and flexible scalability, we applied the fLTRG 
algorithm to a spin- 1/2 frustrated bilayer honeycomb Heisen- 
berg model to which the QMC is not directly accessible, 
and disclosed a quantum phase transition (QPT) from an or- 
dered antiferromagnetic phase to a gapless quantum spin liq- 
uid (QSL) in the ground state. The thermodynamic properties 
of this frustrated spin system are also calculated. Our results 
manifest that the fLTRG would be very promising to tackle 
the intractable correlated quantum many-body systems in two 
and higher dimensions. In what follows, we shall describe the 
basic procedure of the fLTRG algorithm with a quantum spin 
system as an example on a honeycomb lattice. 

Initialization. — Suppose that the Hamiltonian of the sys- 
tem can be written as H = J^uHij, where H,j is a local 
Hamiltonian of pairs of spins. The partition function Z is the 
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FIG. 1: (Color online) (a) The local evolution operator U'L, is de- 
composed via an SVD into two gates, Gfj s and Gf,„ s , each of which 
has two physical bonds (i, i' and _/,/, black) and one geometrical 
bond (s, blue); (b) Contract the shared physical bonds among G L and 
G R to get tensors T L and T R ; (c) A TPDO with inverse temperature 
t. Note that the singular value vectors A 1 " '" on each geometrical 
bond is not indicated in (b) and (c) for concise. 



trace of the density matrix p = exp( —f}H) with ft = l/T the 
inverse temperature and £g = 1. By means of the Trotter- 
Suzuki decomposition, the density matrix can be written as 
p a [exp(-TZijH u )] K+ \ where /3 = (K + l)r, and r is 
the infinitesimal imaginary time slice. Define a local evolu- 
tion operator C/y = exp(-T#y). Then, the density operator 
can be represented as p [fl u Uij] K+l = nS' Ilij #y> 
where f is the Trotter index. In this way, the density matrix 
p is transformed into a TN. By making a singular value de- 
composition (SVD) on Uf., = (ij\Uij\i'f) where stands 
for the direct product basis of spins at site i and j, we have 
[/!/, = Tj,G L ., A { lG R ., , where A is the singular value vector, 
and G L and G R are two local evolution tensors, each of which 
has two physical bonds (i, i' and 7,/, respectively) and one 
geometrical bond (s). For a honeycomb lattice, this step is de- 
picted in Figs. [T|(a) and (b). Next, by contracting the shared 
bonds among G L and G R [Fig. [U(b)], we get 

1 U,xyz ~ Zu jk,y U kl.z' 1 il,xyz ~ /j Cr y,* Cr jky^kW ( 1 ^ 

jk jk 

where x, y and z are three inequivalent bonds on a honeycomb 
lattice [Fig. [T](c)]. The density operator p at an inverse tem- 
perature t has the form of 
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in which Trg is the trace over all contracted geometrical 
bonds, and A 1 , A", A IU are three inequivalent singular value 
vectors with the initial value A . This gives a tensor prod- 
uct density operator (TPDO), which is a direct extension of 
the matrix product density operator [2 1 ] and the tensor prod- 
uct states. In fact, the TPDO is the infinite product of two 
inequivalent tensors T L and T R for two sublattices (denoted 
as <SX fl and SJLb) of the honeycomb lattice as well as A 1 , A 11 
and A 111 for three inequivalent bonds [Fig. Q](c)]. Because of 
the structure of the present lattice and the forms of interac- 
tions, only two inequivalent tensors are adequate here, which 
is independent of any specific states. For a Kagome lattice, at 
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FIG. 2: (Color online) (a) A TPDO on bond x; (b) acting G L and G R 
on (a); (c) contracting the non-orthogonal transformation matrices 
V L and V R ; (d) obtaining new tensors T L , T R and matrix M; (e) 
decomposing M by employing the SVD and contracting the top D c 
left and right singular vectors T L and T R , respectively, to obtain the 
truncated T L and T R \ and (f) keeping the D c largest singular values 
as new A 1 . 



least three inequivalent tensors are needed. We now present 
the fLTRG process on bond x [Fig. |2](a)] as an example. 

Evolution. — By acting G L and G R in pairs on the TPDO to 
evolve along the imaginary time direction, we have 
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as shown in Fig. [2] (b). The dimension of the newly gained 
bond (xx') is enlarged. We denote (xx') as an index a. Then, 
the corresponding singular value vector A 1 is a direct product 
of A and A 1 : A' a = A° X A' X ,. To obtain the optimal approxima- 
tion of the TPDO with the truncation of the enlarged bond, 
we propose a modified HOOI (mHOOI) algorithm. The origi- 
nal HOOI takes the interactions among each mode of a tensor 
into account by iterating the orthogonal transformation on ev- 
ery bond, implying that the effect of other bonds (or modes) is 
thus considered when truncating one bond. For our purpose, 
we suggest the mHOOI algorithm by considering the interac- 
tions of not only each bond but also each tensor. 

mHOOI.— Define the reduced matrix M L (M R ) of f L (f R ) 
on bond a by 

K = £ f iAzW H , *J = 2*^5*^^.(4) 
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Making an eigenvalue decomposition on M L and M R , we have 

M L ap = A^A^, M R p = £ A^r*A^, (5) 

x x 

where the matrix A L (A R ) is formed by the eigenvectors of 
M L (M R ), and F L (T R ) contains the corresponding eigenval- 
ues. The non-orthogonal transformation matrix V L (V s ) can 
be obtained by 



V L =A L A /F V R =A S ^[r* 

"X <*X y X' a X "X V X 



(6) 
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Acting V L and V R to A l a and their inverses to T L and T R , re- 
spectively, as shown in Figs. |2](c)-(d), one has 
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It corresponds to inserting two unit matrices [/ = (V L Y 1 ■ V L — 
V R ■ (V R y l , as shown in Fig. [2] (c)] and changes nothing for 
the TPDO. The intermediate matrix At can be decomposed 
through SVD as 
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where ~A K is the singular value vector arranged in a descending 
order, P (Q) is formed by the left (right) singular vectors of At. 
Now, we keep the D c largest singular values as new singular 
value vector A 1 of bond x, and normalize A 1 by dividing the 

renormalization factor = ^^^(/H) 2 with n the step of 
evolution. Meanwhile, acting the top D c singular vectors in P 
and Q on T L and T R , respectively, we get new tensors with 
truncated bond x [Figs. |2](e)-(f)] 
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Then we renew A" , A 1 " and A 1 in turn without truncating their 
dimensions by making the iteration procedure several times 
(e.g. five times in our case) according to the operations de- 
scribed in Figs. |2](c)-(f) until reaching a convergence. 

fLTRG step. — The evolution and mHOOI processes give a 
complete fLTRG step on bond x. Doing this step on x, y and 
z bonds in one turn corresponds to that the TPDO is evolved 
with an imaginary time r. After doing the Kth-Xxmi, the in- 
verse temperature for the TPDO reaches — (K + 1)t. Con- 
sequently, the density operator p is obtained by Eq. (f2]). 

It should be remarked that in the above mHOOI procedure, 
we first make the truncation on bond x and then do the itera- 
tion over three bonds so that the interactions among bonds and 
tensors are well taken into account. Certainly, one may also it- 
erate first and then truncate the enlarged bond, which gives al- 
most the same result according to our calculations. However, 
doing the truncation first is obviously more efficient. More- 
over, for the present case with an infinite size, we have only 
three inequivalent bonds on which the iteration goes. In prin- 
ciple, such an mHOOI may also be applied to the finite-size 
systems by sweeping over all inequivalent bonds to achieve 
the optimal approximation. 

Free energy. — Partition function Z can be obtained by trac- 
ing all physical and geometrical bonds. Tracing all physical 
bonds of the TPDO, we get a 2D classical TN (CTN). The 
free energy per site / = - lim^K, In Z(J3)/ (N/3) with N the 




FIG. 3: (Color online) The inverse temperature fi dependence of (a) 
the energy per site E for different 6 with h s = and (b) the staggered 
magnetization per site m s for different h s with 5 = 0.5 in a spin- 
i /2 anisotropic Heisenberg antiferromagnet on a honeycomb lattice, 
where D c = 22 and r = 0.005. The QMC results are also included 
for comparison. 



number of lattice sites, is comprised of two parts, the renor- 
malization factors r^ and the contributed factor per site ru 
obtained through the contraction of the CTN: 
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(12) 
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The thermodynamical quantities including energy, magnetiza- 
tion, susceptibility and specific heat of the 2D quantum sys- 
tems can thus be obtained. 

What is more, the ground state properties can also be stud- 
ied with the fLTRG algorithm. When one takes K — > oo and 
t — > 0, the renormalization factors of each fTRG step con- 
verge to 1 . The ground state energy per site eg nas a simple 
form of 
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Spin- 1/2 Heisenberg antiferromagnet on honeycomb lattice. 
— To test the efficiency and accuracy of the fLTRG algorithm, 
we employ the spin- 1/2 anisotropic Heisenberg antiferromag- 
net on honeycomb lattice in a staggered magnetic field h s , and 



compare the fLTRG results with the QMC calculations 122 1. 
The local Hamiltonian of nearest-neighbor spins reads 



■S))h,/3, 



(14) 



where S x , § y and § z are the x-, y- and z-component of spin 
operator on the zth site, respectively, and 6 measures the 
anisotropy of spin couplings. The energy per site can be cal- 
culated by E = -d(J3f)/d/3, and the staggered magnetization 
per site is obtained by m s = df/dh s . Fig. [3] gives E and m s 
as functions of the inverse temperature /3 for different 6 with 
h s = and different h s with 6 = 0.5, respectively. It can be 
seen that our fLTRG results are in nice agreement with those 
of QMC calculations, showing that the fLTRG algorithm is 
feasible, efficient and accurate. 

The specific heat as a function of (3 is calculated by C — 
-p 2 dE/d/3, as shown in Fig. |4]for S = 0.5. A divergent peak 
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FIG. 4: (Color online) The inverse temperature f} dependence of the 
specific heat for 6 = 0.5 and h s = in a spin- 1/2 anisotropic Heisen- 
berg antiferromagnet on a honeycomb lattice. The QMC result is 
included for a comparison. Inset gives the 6 dependence of the criti- 
cal temperature T c . Here D c = 32. 



at a critical temperature T c is observed, which indicates that a 
phase transition occurs between a paramagnetic phase and an 
antiferromagnetic phase at T c . It is also well consistent with 
the QMC result, showing again the efficiency and accuracy 
of the fLTRG method. In the inset of Fig. [4] T c as a func- 
tion of 5 is given, indicating that T c declines almost linearly 
with increasing 6. It should be pointed out that as 5 — > 1, the 
divergent peak of the specific heat becomes gradually round 
owing to the increase of quantum fluctuations, and the phase 
transition no longer exists at 8 — 1, being consistent with the 
Mermin-Wagner theorem. 

Spin- 1/2 frustrated bilayer honeycomb Heisenberg model. 
— To show the power of the fLTRG algorithm, we now apply 
it to investigate the spin- 1/2 anisotropic Heisenberg antifer- 
romagnetic model on a bilayer honeycomb lattice with alter- 
nating antiferromagnetic and ferromagnetic interlayer interac- 
tions J a and J/, [the inset of Fig. |5](a)]. It is a quasi-2D quan- 
tum frustrated spin system to which the QMC is hardly acces- 
sible owing to the negative sign problem. The local Hamilto- 
nian of this model is defined as 



Ha = H 



in 



■ &ff + (Hf ] + #f)/3, 



(15) 



where tiff = J y [S y (S*S X j + + S^Sj], with the layer 

index y — 1 and 2, and i, j the nearest neighbor sites within 
the single layer; tif M = J a , b [6 a , b (§ fS * + S y .$*) + is 
the interlayer couplings. When J a > and Ji, < 0, it gives 
rise to the spin frustration. Without losing generality, we shall 
take Ji = J 2 = J > 0, J a > 0, J b < 0, J a = -Jb = /', 
and 5\,2 = 6 a ,b = 8 = 0.5. As the frustration exists, this 
model would be expected in proper circumstances to have a 
QSL ground state that is currently under an active debate 123 - 

m 

Figure [5] (a) shows the sublattice magnetization per site m z 
as a function of the coupling ration /' /J at zero temperature. It 
can be seen that there exists a quantum critical point (J' /J) c = 
2.6, at which a QPT occurs. When J' /J < 2.6, m z is nonzero 




FIG. 5: (Color online) (a) The sublattice magnetization per site m ; 
as a function of J' /J for a spin-1/2 frustrated bilayer honeycomb 
anisotropic Heisenberg model (inset) with 6 = 0.5 at zero tempera- 
ture, where a quantum critical point (J'/J) c = 2.60(2) is identified, 
(b) The inverse temperature fS dependence of the specific heat for 
different J' /J with h s = 0. Here D c = 32. 
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FIG. 6: (Color online) (a) Log C versus /? of the spin-1/2 frustrated 
bilayer honeycomb Heisenberg model for various J'/J with h s = 0; 
(b) The gap A versus J'/J, where the extrapolation shows that at 
(J' IJ) C = 2.60(2), the gap vanishes, indicating a QPT to a gapless 
QSL state. Here D c = 35. 



and decreases slowly with the increase of J'/J, showing that 
in this regime the system is in an antiferromagnetic ordered 
state; around J'/J 2.6, m z drops sharply; and it goes to zero 
for J'/J > 2.6, while the magnetization in the x - y plane 
is found about 10~ 3 in the whole region, suggesting that the 
system enters into a disordered state. This disordered state is 
nothing but a gapless QSL state (see below). The reason is 
that, for J'/J > 2.6, the frustration becomes stronger If26ll . 
which strongly suppresses the magnetic long-range ordering, 
giving rise to a QSL state. 

The temperature dependence of the specific heat (C) of this 
frustrated spin system is given for different J'/J in Fig. 0(b). 
When J'/J < (J' /J) c = 2.6, the specific heat displays a diver- 
gent peak at a critical temperature for a given J'/J, showing 
a second-order phase transition between a paramagnetic state 
and an Ising-type ordered state. It appears that the critical 
temperature depends weakly on J'/J. For J'/J > (J' /J) c , 
the specific heat shows a round peak, and no phase transition 
happens, which is consistent with the observation that in this 
regime the ground state of system is in a gapless QSL state. 

Finally, we observe that, when J'/J < 2.6, this frustrated 
bilayer spin model is in an Ising-type ordered state with a gap. 
It is evidenced by the low-temperature behavior of the specific 
heat that decays exponentially with the inverse temperature fi 
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in a form of C ~ exp(-A/T). The gap A can be determined 
by using a linear fitting between log C and {S for different J' I J 
with h s = 0, as presented in Fig. |6](a), showing a perfect linear 
J' I J dependence of the gap [Fig. |6](b)]. By extrapolation, one 
may observe that A vanishes at (/' /J) c = 2.60(2), confirming 
the QPT from a gapped Ising-type ordered state to a gapless 
QSL state when the frustration effect becomes stronger. 

In conclusion, by extending the Tucker decomposition to 
a tensor network, we propose a novel algorithm coin as the 
fLTRG, and examine its efficiency and accuracy by employ- 
ing a spin- 1/2 Heisenberg antiferromagnet on honeycomb lat- 
tice. The fLTRG results are well in agreement with the QMC 
calculations. To show the power of the fLTRG algorithm, it 
is applied to a spin- 1/2 frustrated bilayer honeycomb Heisen- 
berg model with alternating interlayer couplings, and a quan- 
tum phase transition is disclosed, where a gapless quantum 
spin liquid phase is identified. The present fLTRG algorithm 
could be straightforwardly extended to other quantum lattice 
systems. 
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